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Abstract 

Frequently, during the process of solving a mathematical model numerically, we end 
up with a need to operate on a vector v by an operator which can be expressed as /(A) 
while A is N x N matrix (ex: exp(A), sin(A), A -1 ). Except for very simple matrices, it is 
impractical to construct the matrix /(A) explicitly. Usually an approximation to it is used. 
In the present research, we develop an algorithm which uses a polynomial approximation 
to f(A). It is reduced to a problem of approximating f(z) by a polynomial in z while 
z belongs to the domain D in the complex plane which includes all the eigenvalues of 
A. This problem of approximation is approached by interpolating the function f(z) in a 
certain set of points which is known to have some maximal properties. The approximation 
thus achieved is “almost best.” Implementing the algorithm to some practical problems is 
described. 

Since a solution to a linear system Ax = b is x = A -1 6, an iterative solution to it can 
be regarded as a polynomial approximation to /(A) = A -1 . Implementing our algorithm 
in this case is described too. 
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I. Introduction. 


Let A be a N x N matrix whose elements belong to C, and f{z) is a function such 

that 


The matrix f(A) 

f(z) : C . 

can be defined in the following way: If 

(1.1) 


OO 

f( z ) = ^“iiz-zoY 

\z — Zo\ < r 

(1.2) 

then 

1=0 

OO 

f(A) = 'E*i(A-i oiy 

t=0 

|A k - z 0 \<r 

(1.3) 


where A* is an eigenvalue of A. In (1.3), f(A) is expressed as an infinite polynomial. It can 
be shown [Gant 59] that / (A) as defined above can be written also as a finite polynomial 
of degree < N — 1 as follows: 

f(A) = £ + — + (1.4) 

k= 1 

where 

Afc — an eigenvalue of A 

jk — multiplicity of A* (1.5) 

Hkj{A) — polynomial in A of degree < N — 1 . 

In many practical applications, we would like to compute a vector u> which results 
from operating with the matrix f[A) on a vector v 

w = [/(A)]v . (1.6) 


Using expression (1.4) for this purpose has two major disadvantages: 

(a) The exact knowledge of the eigenvalues is required. 

(b) / (A) is expressed as a polynomial of degree < N— 1 which can be large; thus it results 
in an highly time consuming operator. 

In this paper, we would like to present an algorithm which approximates the matrix f(A) 
by a polynomial P m (A), where m << N. Since (1.4), we have 

E m (A) = f(A ) - P m (A) = £ [£(A*)H M (A) + ... + £<*-» {X k )H kii (A)] (1.7) 

/c=l 


I 
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while 

E(z) = f(z) - P m (z) . (1.8) 

Hence, our problem is reduced to a problem of approximating the function f(z) by a 
polynomial P m (z ) while z belongs to a domain in C which includes all the eigenvalues of 
A. It is obvious from (1.7) that in the case where > 1, Pm k ~ ^ has to approximate 
yOfc-i) a t the point A*. 

In Section 2, we describe briefly how to approach this problem of approximation. Using 
the results of this section we can construct a computational algorithm, once a suitable 
conformal mapping function is known. In very few cases, it is possible to get this function 
analytically. Usually we have to resort to a numerical method. A lot has been done in this 
area of conformal mapping and there are suitable routines. In our work, we have used a 
set of routines written by N. Trefethen based on his paper [Tref80], which are very efficient 
and reliable. The interested user can get the routines from the author upon request. The 
algorithm which results from making use of the approximating polynomial is presented in 
Section 3. In Section 4, we deal with the rate of convergence of the suggested method. Few 
examples of mathematical models for which the new algorithm can be implemented are 
given in Section 5. The numerical solution of a system of linear equations Ax = b where 
A is a general nonsymmetric matrix can be regarded as a particular case of our problem 
where the function which has to be approximated is f(z) = 1 jz. Section 6 is devoted to this 
subject. Using tools from the theory of approximation in the complex plane for inverting 
nonsymmetric matrices has been studied already by other researchers. A few of them 
are: T. Manteuffel [Mant77], [Mant78], D. Smolarski and P. Saylor [SmSa85], Gutknecht 
[Gutk86], Y. Saad [Saad87],and others. The algorithm discussed in our paper gains its 
simplicity from the fact that it is based on the powerful tool of interpolation. Thus, it can 
be implemented in a straightforward way, once the conformal mapping function is known. 
An important conclusion of the analysis is that for the general nonsymmetric case, the 
significant factor which governs the rate of convergence of most of the iterative algorithms 
is not the well known condition number ||A|| • ||A -1 ||. It is shown that the relevant factor 
is p/R where p is a parameter which measures the size of the domain and R is related to 
the orientation of the domain with regard to the singular point of the function \ (z = 0) . 
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Based on this conclusion it is shown that some iterative algorithms, like standard SOR, 
are based on an incomplete optimization process. We conclude this paper in Section 7 by 
giving some numerical examples. 

2. Polynomial Approximation in the Complex Plane 

Let D be a bounded closed continuum in C such that the complement of D is simply 
connected in the extended plane and contains the point at oo. Define now 

A(D) - the space of functions which are analytic at every point of D. 

7r m - space of complex polynomials of degree < m. 

Then it is well known [SmLe68] that for every / 6 A(D) there exists 6 7 r m such 

that 

||/-PAIIoo<||/-Pm||oo VP m £ 7T m . (2.1) 

From algorithmic point of view, it is relatively complicated to find P In many cases, it 
is quite simple to find a polynomial approximation which is “almost” as good as P^. It is 
found by projecting A(D) on 7 r m . If S m is a projection operator 



Sm • > 7T m , 

(2.2) 

then 


(2.3) 

thus 

11/ - S m (/)|| < (1 + ||S m ||)||/ - P^|| . 

(2.4) 

If ||5 m || is reasonably small, regarding S m [f) as “almost” as good as P^ 
example, if D is the unit disc and 

is justified. For 


s.m-t'"?-': 

}=o J - 

(2.5) 

then [GeMa75] 

< — log(m) + 1 + 0(1) . 
7 r* 

(2.6) 
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If D = [—1, 1] and 

m 

s.(/) = E** r ‘( 2 ) < 2 - 7 ) 

k= 0 

while 

Tk(z) = cos(k(cos~ 1 z)) (Chebyshev polynomials) , (2.8) 


then the bound on ||5 m || is exactly the same as in the previous case (2.6). A generalization 
to an arbitrary domain D (as defined in the beginning of this section) is given by a Faber 
polymomials expansion. Let 4>(z) be a conformal mapping from z to w which maps the 
complement of D to the complement of a disc of radius p such that 


lim M 

z-*oo z 


= 1 . 


(2.9) 


p is the logarithmic capacity or the transfinite diameter of D [SmLe68]. 
If the Laurent expansion at oo of [<j>(z)] m is 


[<f>(z)] m = z m + c m -iz m 1 + • • • C\Z + cq + C—\Z 1 + ' -- , (2-10) 


then Faber polynomial of degree m, F m (z ), which corresponds to the domain D is the 
polynomial part of (2.10): 


fmW = z m + c m -i z m 1 H 1- Co ■ 


(2.11) 


We have 


Cj 2 7T* L, =r 


i m: 

ZJ + 1 


-dz 


( 2 . 12 ) 


while R is chosen sufficiently large so that D is contained in the interior of the region 
bounded by the circle \z\ = R [SmLe68]. Given / € A(D ) then 


/(*) = ■ 


(2.13) 


k—0 


dk 


i f fjjj M l 

27 ri J | w |_ r w k+1 


dw 


The coefficients a& are 



ORIGINAL r,-.. 

OF POOR QUALITY 


where ip(u>) is the inverse of <j>(z) and r > p is sufficiently small that / can be extended 
analytically to I r . ( I r is the closed Jordan region with boundary T r where T r is the image 
under xjj of the circle |w| = r.) In [Ella83] it was shown how the cy’s and s can be 
computed efficiently using F.F.T. Based on (2.13), the matrix f{A) can be approximated 
by P m {A) where 

m 

P m (A)=£a*F*(A) • (2.14) 

k= 0 

When D is an ellipse with axes parallel to the real and imaginary axes, then Fk is the 
translated and scaled Chebyshev polynomial T* [Mark77]; thus it satisfies a simple three 
term recurrecnce relation. This recurrence relation enables us to compute 


u = 


m 


Y. akFk{A) 


k=0 


(2.15) 


efficiently. A detailed description of these cases, approximating the operator e A , can 
be found in [Tale85], [Tale86], [KoTa86], [TaKo84]. For more complicated domains, Fk 
satisfies a k terms recurrence relation [Ella83]; thus from memory point of view it is 
impractical to use an algorithm based on (2.15). 

This major drawback can be overcome by using a different approach to the approxi- 
mating problem. It uses interpolation as the projection operator S m . Using this tactic we 
are confronted with the following question: Which are the interpolating points in D, such 
that the interpolating polynomial will be “almost” as good as It is known [SmLe68] 
that if D is the unit disc, then two possible sets of points are: 

(a) Zj = 0 j = 1, ... ,m (zeroes of z m ) 

(b) The m roots of unity. 

In a similar way, for a general domain D, a two possible sets of points are [SmLe68]: 


(a) the m zeroes of F m (z) (2.14a) 

(b) Zj = ip(wj)j — 1, . . . , m, while uij are the m roots (2.14b) 

of the equation w m = p. 

In the first case, it can be shown [GeMa75j that 


ORIGINAL PASS m 

0E POOR QUALITY 

while in the second case 

H-S'mll < ^log(m) + 1 + 0(1) . (2.16) 

Since, interpolating at Zj = J = 0 , . .. ,m does not involve computing Faber poly- 

nomials, it is the simplest and most efficient approach to this problem of approximation. 
The interpolating polynomial in Newton form is: 

P m {z) = a 0 + ai(z - z 0 ) + a 2 {z - z 0 )(z - zi) H 1- a m (z -z 0 )---(z- z m - 1 ) (2.17) 

where a* is the divided difference of order k 

ak - f[zo,...,z k ] A; = 0, . . . ,m . (2.18) 

When p ( the logorithmic capacity of D) is large, will be very small; thus it is preferable 
to make the following change of variables 

z = z/p . (2-19) 

Hence 

f{z) = f(z) = f(p • z) (2.20) 

and 

P m (l) ^ f(z) (2.21) 

where 

Pm{z) = b 0 + bi(z - z 0 ) H h b m {z - z 0 ) ■ •• (z - £ m _i) (2.22a) 

bk = /[^ 0) • • • > k = 0, . . . ,m . (2.226) 

The only difficulty in finding P m (z) is to get the function V’( w )- For simple domains, 
it is possible to find this function analytically. For more complicated domains one has to 
resort to a numerical approach. 

When the domain D is a polygon, the mapping function is Schwartz-Cristoffel trans- 
formation. In [Tref80], a very reliable and efficient algorithm for mapping the interior of 
the unit disc to the interior of the polygon is described. Since, in our case, the mapping 
of the exteriors is needed, the routines should be modified accordingly. 
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Separated, domains: In some very important physical phenomenon, we have a situation 
where the eignevalues of the matrix A are clustered in two or more separated domains which 
can be far apart (for example: stiff O.D.E.’s problems of parazitic waves in meteorological 
sciences, spectral methods for nonperiodic boundary values problems, etc.). Hence, the 
domain D is a union of k subdomains: 

D = uLiA- • 

In this case, the complement of D is not simply connected any more but just connected. 
The basic theory regarding the simply connected case extends to the more general one. A 
detailed analysis of this problem is carried out in our next paper. It is shown there that 
the interpolating points can be taken as a union of sets of points achieved by considering 
each domain separately. In Section 7, we bring some numerical examples of this case. 


3. The Algorithms 

Based on (2.22), we will approximate the operator f{A) by P m (A) while 
P m (A) =b 0 I + b 1 {A- z 0 I) + b 2 (A — z 0 I){A - zj)+ 

h bm{d — zol) • • • (A — Z m - 1 1) . 

A = (1/p) A 

Operating with P m (A) on a vector v is carried out by the following algorithm 
Algorithm 1: 

u = v 
<*> = bov 

For i = 1,. ..,m do 
u = (A — Zi-\I)u 
xv = w + b{U 

The output of Algorithm 1 is the vector w which satisfies 

U> = P m {A)v . 


(3.1) 


(3.2) 


(3.3) 


Roundoff errors of Algorithm 1 depend strongly on the arrangement of the interpolating 
points. 
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In Appendix A we describe how to arrange the points, taking into account this phe- 
nomenon of roundoff errors. 

In many practical problems, we have the following situation: 

1) A is a real matrix (3.4a) 

2 ) /(*) = /(*)• ( 3 - 4b ) 

In this case, it is possible to design an algorithm similar to (3.2) which will be carried 
out in real arithmetic even so Z{ and 6,- are complex numbers. This result is based on the 
following two theorems: 

Theorem 1: Let z 0 ,zi,z 2 ,z 2 , ... ,z k ,Zk be 2k interpolating points where z 0 and z x are 
real numbers. If P 2k -i{z), 

Pik-i {z) — a o + a i z + ' • • a 2 k-iZ 2k 1 (3.5) 


is the interpolating polynomial of a function f(z) which satisfies (S. 4b) then 
ao,ai,.. .,a 2 fc-i are real numbers. 

Proof. The Lagrange formula of P2k-i(z) is 


k k 

P2k-i(z) = Y £j{z)f{zj) + Y t)(z)f(zj) (3.6) 



04 

11 

o 

II 


where 

<;(*) = §TJ-\ — - — j = ° k 

Qjk(zj) z ~ z i 

(3.7) 


«;(*) = i = 2,...,k 

3 Q2k( z i) z ~ z i 

(3.8) 


k k 

Q2k{z ) = ~ *»') II ( 2 “ *»') * 

(3.9) 


i=0 t=2 


Since (3.9) we have 

k k 

Qu ( z j) = n & - n (*> ~ ^ y = • ■ • » * 

(&$) 


(3.10) 
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Hence 


Therefore, 


k k 

Q2ki z i) ~ Yl( Z J ~~ 2 ») JJ { z j ~~ z i ) J = 2, . . . yk . 

(£) 


^2fc( z y) — II ( z y z ») Il( z y — z ») — 

( 535 ) 

*: k 

— JJ ( z y — z i ) JJ( z y — z ») = Q2k( z j) 

( 535 ) 

k k k k 


i=0 i=2 

= Q2fc( 2 ) • 


/ / \ ^2fc( 2 ) l _ f ^y ( z ) 

<%*(*,)»-*, U*(*) 

T+TZT Q2fc(^) 1 


Using (3.4b) and (3.14) in (3.6) we get 

k 


P 2 fc-l( 2 ) = e*{z)f( Zj ) + ^2 ^( z )/( z j) = ^ > 2k—l ( 2 ) 
y =2 y=o 


(3.11) 


(3.12) 


II 

sa 

- z i) ]3(* “ *»') = 

i=0 

(3.13) 

1bT1*r 

* 'r> 

j = o, fc 
2 < j < 1 

(3.14a) 

= ^y( z ) 

2 <j<k. 

(3.146) 


(3.15) 


and the proof is concluded. 

Theoretically, it is possible to write an algorithm based on (3.5). It means to approx- 
imate f(A) by Pm(A) where 


P m (A ) = aol + ai A + • • • + a m A m . 


(3.16) 


This algorithm cannot be used for practical problem since huge roundoff errors result. We 
would like to stick to the Newton-form which is much more robust from the roundoff errors 
point of view. For this purpose, we state and prove the following theorem: 

Theorem 2: Let P 2 fc-i( 2 ) be the interpolating polynomial as defined in Theorem 1, written 
in Newton-form 

k- 1 

P 2 k-i{z) = b 0 + b 1 {z- z 0 ) + {z- z 0 ){z - z^ Si{z)Ri{z) 

i=l 


(3.17) 
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Si(z) — 6 2 j + t>2i+\(z — Zi+ 1) 


Ri{z) = n ( Z -Zj){z-Zj) » = 2,...,fc-l 


(3.186) 


-Ri(*) = 1 • 


(3.18c) 


P 2 *-i(*) = 6 0 + 6 x (z - z 0 ) + (* - «o)(« - *i) II S *i z ) R i( z ) ( 3 - 19 ) 


S*(*) = b% + 6f t+1 (z - afti) 


(3.20a) 


6& = Real (6 2t ) 


(3.206) 


^ 21+1 — Real (6 2t -+i) — ^ 2 *-i-i 


(3.20c) 


*/+ i = Real (z»+i). 


Proof. We have 


iWl(*) = -P2fc-l(«) + (^ - 2 0 )(2 - *l)Sfc(z)R*(z). 


Define 


7r - The set of all polynomials with real coefficients. 


Then, by theorem 1 


Rk{z) € 7r by definition. 


-P 2 lc+l(z)> P 2 k-l{z) <£ 7T • 


Since 


5k(«) = [P 2 fc+l(*) - P2k-x{z)\l{z - Z 0 )(z - Z\)Rk{z) 


s k e 7T 


Thus 


^2k+l — l>2k+l 
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&2k — hk+lZk+l — &2k — ^2k+l z k + l (3.24) 

and the proof is concluded. 

Based on (3.19) , the vector uj 

u=[P 2 k-i{A) }v (3.25) 

is computed by the following algorithm: 

Algorithm 2. 

w = [b 0 I + b x (A — ZqI)]v (3.26a) 

r = (A — zoI)(A — zil)v. (3.2 6b) 

For t = 1, , . . . , fc — 1 do: 

r = (A — zf +l I)r (3.26c) 

w = w + b%r + 6f <+1 f (3.26 d) 

r = (A - zf+ x I)r + ( z/ +1 ) 2 r {zj +1 = I m (z i+1 )). (3.26c) 


This algorithm requires three vectors. It is possible to save one vector by using an 
algorithm based on (3.16). As mentioned previously, this algorithm has the disadvantage 
of sensitivity to roundoff errors. Thus, it can be used only with low degree polynomials. 
Another possible way of saving one vector and which is much more robust from the roundoff 
errors point of view is to apply (3.19) through calculating the roots of F > 2k-i( 2 )- Since 
P 2 k-\{z) is a polynomial with real coefficients, we have the following set of roots 

Al, A 2 , • • • , A r , /li, Ml 5 • • ■ ) Ma> Ms (3.27) 


where 

r + 2s = 2k-l] A, are real . (3.28) 


Hence 


P 2 k-i(A) = a - A il) f[(A 2 - 2 Rem A + \m\ 2 I) 


1 = 1 


1=1 


a — 


i , 2fc-l(0) 


r s 


p = (- i ) r ii A *n n 2 • 

i= 1 *=1 


(3.29) 


(3.30) 
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Operating with (3.29) requires only two vectors. This approach is limited mainly by the 
sensitivity of the algorithm which finds the roots of P^k-iiz). 

4. Rate of Convergence. 

We have the following definitions: 

Definition (4.1): Let T# be the image under of the circle |u;| = R(R > p ) and Ir is 
the closed Jordan region whose boundary is T.R. If f(z) is single valued and analytic on 
Ir, then the sequence of polynomials P m (z) is said to converge to f(z) on Ir maximally if 

\f{z)-P m (z)\<C(p/R) m z € Ir (4.1) 

where C depends on p/R but not on m or z. 

Definition (4.2): The set of interpolating points Zj = ip(u>j) is said to be uniformly 
distributed oxiTd (the boundary of D) if are equally distributed on the circle |w| = p. 
Using these two definitions, we can quote the following [Wals56]: 

Theorem: Let D be a closed Jordan region. Let the points (3 lie on the boundary 
Tp of D. A necessary and sufficient condition that the sequence of polynomials P m [z) 
of degree m found by interpolation to a function f(z) analytic on D in the points (3^ 
converges uniformly to f{z) on D is that the points (3 be uniformly distributed on r^. 
If this condition is satisfied, the sequence P m (z) converges maximally to f(z) on D. 

Given a domain D and a function f(z), p and R are defined explicitly and we have 

pjR — the asymptotic rate of convergence. 


If f(z) is an entire function, (4.1) is satisfied for arbitrary R. Using Theorem (1), 
we can expect that the algorithm described in Section 3 will converge very rapidly for 
the approximation of the operator exp (A). On the other end, for the operator A~ x , the 
rate of convergence will depend strongly on the size of D (the parameter p) and on its 
orientation with regard to the singular point z — 0 (the parameter R = |</>( 0 )|). Since the 
set of interpolating points Zj = ip (pe t] ^\ j = 0, ...,m — 1 depend on m, given a 
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desired accuracy e, one has to decide a priory on the degree of the polynomial. Deciding 
on m can be done in the following ways: 

1. Getting the parameters p and R (analytically or numerically) and choosing m such 
that 

(/ P/R) m = e . (4.2) 

2. Computing the error numerically on a set of check points on the boundary for different 
m’s, and choosing m which will satisfy the desired accuracy. 

Using 1 or 2 gives us information on | f(z) — P m (z)|oo- Substituting it for ||/(A) — 
imU)||oo is relatively accurate when A is a normal matrix, since in this case we have 

\\f(A) - P m (A)|| O0 < \f(z) - P m ( 2 )| 00 ||T||- 1 ||T|| 

= I f( z ) ~ Pm {z) |oo 

where T and T~ l are the unitary matrices which diagnolize A. When A is “far” from 
normality, m should be increased by an amount which depends on ||T|| • ||T -1 ||. 

5. Applications. 

Frequently, while solving a system of O.D.E.’s or P.D.E.’s, we end up with the following 
set of differential equations 

~u N -a K u r] = s tl (x,t) (5.i) 

£W°) = K 

where Un and S^r are vectors of dimension N and G x is a N x N matrix. Expanding 
Sn{%> t) as 

k 

S N {x,t) =^2 a j{t)S 3 N {x) 

y=i 

enables us to write the formal solution of (5.1) as 

k 

U N (t) = exp(G„t)V% + 

/= 1 


(5.2) 
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where 

fj(G N t) = I exp (GNT)ctj(t - r)d .t 1 < j < k (5.3) 

Jo 

Uffit) can be approximated by the algorithm discussed in Section 3 where the functions 
which have to be interpolated are 


f 0 (z) = exp(*z) 


(5.4a) 


fj(z) = f exp (zT)aj(t — r)dr 1 < j < k (5. 

Jo 

In the case where (5.1) is originated from a system of hyperbolic P.D.E.’s, the domain D 
is on the imaginary axis (in the constant coefficients case) or close to it (in the variable 
coefficients case). When (5.1) is originated from a set of parabolic equations the domain 
D is on the negative real axis or close to it. In these two cases, the Faber polynomials are 
scaled and translated Chebyshev polynomials. Thus, an efficient algorithm, which makes 
use of the three terms recurrence relations, can be implemented. These two cases are 
described and treated in details in [Tale86], [KoTa86], [TaKo84], [Tale85]. 

In the more general situation, when we have both parabolic and hyperbolic terms in 
the equation, the domain D is more complicated. Let us look at the following simple 1 — D 
equation 

u t = au + bu x + cu xx (a < 0, c > 0) . (5.5) 



If the solution is periodic in space, then a good approximation can be achieved by pseu- 
dospectral Fourier descretization. We get the following semidiscrete representation of (5.5) 


(^iv)t = a Un + b(UN) x + c(Un)xx 


(5.6) 


where 


N 

Un= a * eikX 

k=-N 

is the projected solution. (5.6) can be written as 


(5.7) 


(U N )t = G n U n 


(5.8) 
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while Gn is the finite domain operator 

Gr , = al + b^ + c^ . (5.9) 

The function e tkx is an eigenfunction of G Thus, if Ajt is an eigenvalue of Gn, then 

A k — o- — ck 2 + ibk |&| < N . (5.10) 


Since a < 0 and c > 0, we get that the domain D is the following parabola in the complex 
plane: 

D = {z — x + iy\x — a — cfc 2 ; y = bk |A;| < N} . (5.11) 

In real applications a, b , and c are not constant, but have space dependence. There- 
fore, the domain D will vary accordingly. In section 7 we report on some numerical 
experiments treating this problem. 

In a joint work with Dr. Dan Kosloff and his colleagues, we investigate the implemen- 
tation of the present algorithm to some real geophysical problems which can be represented 
as (5.1). The eigenvalues are scattered close to a T shape domain D 

D = {z\z € D\ U Z >2 ; D\ = [— a,0] ; D-i — [— 


In this case, we have an analytic expression for (Thanks to Nick Trefethen.) The 

conformal mapping tp{oj) which maps the complement of the unit disc on the complement 
of D is: 


Hence 


where 




~(1 + E)(u+ -) + 1 - E 


1/2 


- 4 


E =yf a 2 + 6 2 /6 


\ H u ) = ^( w /p) 

6(1 + E) 


P = 


M = i 


(5.12) 


(5.13) 


(5.14) 


is the logarithmic capacity of D. Numerical results for this domain are presented in section 


7 . 
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A very important area where our method can be implemented is the one of solving general 
nonsymmetric systems of equations. This is topic of the next section. 

6. Linear System. 

The iterative solution to a set of linear equations which can be written in matrix form 

Ax — b (6.1) 

is a well treated problem in the numerical analysis literature. While very efficient algo- 
rithms have been developed for the case when A is a symmetric positive definite matrix, 
the general nonsymmetric case is still a challenging problem. 

In this section, we would describe an iterative algorithm for the solution of (6.1) based 
on the new approach. Since 

x = A~ l b (6.2) 

we can write the numerical approximation X k as 

x k = P k {A)b (6.3) 

where Pk{ z ) is “almost best” approximation to the function z belongs to the domain 
D which includes all the eigenvalues of A. 

In [Mant77], [Mant78], T.A. Manteuffel has described an iterative algorithm based on 
enclosing the eigenvalues in ellipses in the complex plane, thus getting an approximation 
based on scaled and translated Chebyshev polynomials expansion. The algorithm described 
in our paper is more general since it can be implemented to any domain in the complex 
plane. Its advantage will be significant when the discrepancy between the best ellipse as 
defined in [Mant77] and the domain is relatively large. 

A standard approach for solving (6.1) is known as Richardson algorithm. It takes the 
following form: 

x k+1 = x k - a k [Ax k - I) . (6.4) 

E k = x - x k 


If 


(6.5) 
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we get 

E k = I PZ(A))E° 


where 


k—1 


pi (a ) = nu 

j=0 


<*i A ) 


and ay are optimal in the sense that 


( 6 . 6 ) 


(6.7) 


*<&., ft € (6.8) 

(7Tfc is the space of all polynomials of degree k) . We would like to show that an algorithm 
which results from implementating our approach is equivalent to algorithm based on (6.4)- 
(6.8). For this purpose we have 

Theorem 3. Let Tk(z) be a polynomial of degree k defined onC such that Tfc(O) = 1 and 
let Qk- i{z) be the interpolant of the function L a t the roots of Tk then. 


max IP* ( 2 ) I < max IP* ( 2 ) I 

z€D 1 * v n ~ z<=D ' V n 


1-z Qk-i(z) = Tk[z) . (6.9) 

Proof. We have 

Qk — i (^») ~ * = !»•••» k (6.10) 

Z\ 

where Z{ is a root of Tk. Therefore, the polynomial Rk[z) 

R k {z) = T k {z) - (1 - 2 Q k -i{z)) (6.11) 


vanishes at A: + 1 points: 0,zi, . . . , Zk and thus it is identical zero. Hence, 


Tk{z) = 1-2 Q k -i{z) 


and the proof is concluded. 

If x° is the initial guess, we have 


Ax° = b° 


( 6 . 12 ) 


(6.13) 


A(x -x°) = b-b° 


(6.14) 
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Thus 

x = A -1 (6 - b ° ) + x° . 

Approximating A -1 by the interpolant Qfc_i(A) results in 

x fc = Q*_i(A)(6-6°) + x 0 

using (6.5) and (6.14) we get 

E k = \I - AQ k -i(A)}E° . 


(6.15) 


(6.16) 


(6.17) 


Since (6.9), the equivalence is established. 

We have shown that using the algorithm based on interpolating the function j at 
z\ > • • • > Zk will reproduce identical results to Richardson iterations 


x 


,+1 = x J — — b) 


with 


Otj = — 
Z 3 


j = 0, . . . , k — 1 . 


(6.18) 


(6.19) 


Writing (6.18) , (6.19) as an algorithm with real arithmetic takes the following form: 
Algorithm 3 (Richardson). 

(The parameters are ay = l./zy and Zj are defined in Appendix A.) 


4 - „o 


for * = 1, ... , 


- a 0 (Ax° - 6) 

(6.20a) 

- a^Ax 1 — b) 

(6.206) 

Ax 2 * -1 - b 

(6.21a) 

|a,| 2 A]i2* (af = i?e(a,)) 

(6.216) 


Prec onditioning. 

Usually, solving a linear system of equations (6.1) is composed of two stages: 
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1. Modifying the original system to an equivalent one 


Ax = b 


( 6 . 22 ) 


such that (6.22) will converge faster than (6.1). 

2. Solving (6.22). 

In many cases, we have a family of matrices A which depend on a parameter u such 

that 

A„x = b„ (6.23) 

and we would like to choose the optimal w. Based on Section 4, it is obvious that the 
significant factor which determines the convergence is 

A = p/R . (6.24) 


Optimization is achieved when 

|A(w op t)| < |A(w)| Vc o . (6.25) 

Therefore, we conclude that for matrices whose eigenvalues are scattered in the complex 
plane, one should consider the factor p/R rather than the standard condition number 
cond(-A) = ||A|| p -1 ||. It is possible to have two matrices A and B such that 

cond (A) < cond ( B ) (6.26) 

and on the other hand 

Ip/RM > W*m ■ (6.27) 

For example: Let A be a normal matrice whose spectrum r(j4) is 

r(A) = {z\ 1 < \z\ <2; Re z > 0} . 

Let B be a normal matrice whose spectrum is 


r(B) = {z\ \z-3\ < 2} . 
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We have 


On the other hand 


cond(A) = 2 
cond(J3) = 5 . 

[p/R]{A) = 0.8547 
[p/R]{B) = 0.6666- •• . 


(6.28) 


(6.29) 


Since (6.29), Richardson algorithm for B will converge much faster than for A, even so 
cond(B) > cond(A). Based on this discussion, we would like to show also that the standard 
S.O.R. procedure can be considered as an incomplete optimization process. Solving (6.1) 
by S.O.R., we use the following iterative procedure 


x k+1 = T u x k + b u 


(6.30) 


where 

T„ = M~ l Nu 
b w = M~ l b 
M u = D u)L 
Nu = (1 — w)D — u>U . 

U, D, L are the lower, diagonal, and upper parts of the matrice A respectively. The 
optimal w sor is chosen such that 


|r(r„.. r )| = minimal . 


(6.32) 


(r(A) is the spectral radius). 

If x is a solution of (6.30), then 

x — 

(6.33) 

where 

i 

ii 

3 

(6.34) 


Thus, it is evident that rather than using an optimization procedure based on (6.32) 
one should use the more general one, based on (6.25). A well-known example treated in 
the literature for demonstrating the features of S.O.R. is the problem of solving Laplace 
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equation in a rectangle. In [Youn71, p. 205], it is shown that for a certain discretization 
of the rectangle, the optimal u> is 

Cc*3or = 1.25 • 

In this case, all the eigenvalues A{ of T w satisfy 

| A* | = 0.25 . 

Hence, the rate of convergence is 0.25. Optimizing with respect to (6.25) can result with 
a different solution. We do not intend to carry out this optimization process, but to point 
out a different possible parameter w. For this problem, we can choose U 

S7 = 1 + e 0 < e < 0.25 

such that the domain D which includes all the eigenvalues of A „ will be 

D = Bi U B2 


where 


and 


B x = {A 0 } 

i ?2 = {z\ \z — 1| < e} 

0 < Aq < 1 — e . 


Since the complement of D is not a simply connected domain, it is not included in the 
theory discussed in this paper. Implementing our approach to these types of domains will 
be carried out in a future paper. For the time being, let us mention that the basic results 
extend. Thus, since B x is composed of only one point and B 2 is a circle of radius e around 
1, it can be shown that the rate of convergence p/R is 


p/R{A&) = e< 0.25 

and Ao will enter into the constant C (4.1). The set of interpolating points is the union of 
sets of both domains. Since in B x we have only one point, we will get the following points 


A*i * • • • » A^n 
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while fij can be chosen as equally distributed on the circle \z — 1| = e 

Hj = e exp(27rtj/n) j = 1, .. . , n 

or as the n zeros of Faber polynomial of degree m which is ( z — l) n ; thus 

fj,j = 1 j 1 , . . . , n . 

The efficiency of using an algorithm based on U rather then w 30r depends on the accuracy 
needed, since the constant C has been increased by a factor of 1/Ao. 

7. Numerical Results. 

I. Time dependent problem-Parabolic type. 

In this subsection we present numerical results for the following problem 

ut — Gu = S(x,t) (7.1a) 

u(x, 0) = 0 . 

Where 

G = a ^§^ + b ^§^ + C ^ ( 7 * 16 ) 

S ( x , t) = sin (kx) — t { [— k 2 a(x) + c(x)] sin kx + kb(x) cos kx } . (7.1c) 

The exact solution of (7.1) is 


u(x,t) = t sin kx 


(7.2) 


In order to solve (7.1) numerically, we first approximate the space operator G, by the finite 
difference operator Gfj. Assuming periodicity we have 


(G„ «),• = «(*,) "' +1 u ,,+ L - y ~ 1 + i = o, • • • . if - 1 


where N is the number of grid points and 


2Ax 


(7.3) 


Ax — 2ir/N 


(7.4) 
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Xj — j • Ax j = 0, . . . , N . (7.5) 

Hence G ^ is a N x N matrix. The semidescrete representation of (7.1) is 

(«jv)i — Gnun = S + tSpj (7.6 a) 

(u° n ) 3 =0 j = 0,...,N (7.66) 

{S h) j = sin (kxj) j = 0,...,iV (7.6c) 


= \k?a{xj) — c{xj )\ sin kxj — kb{xj) cos kxj j = 0, . . . , N. (7.6 d) 
The formal solution of (7.6) is 

u N {t) = fi(G N t)Sk + f 2 (G N t)S 2 N (7.7) 

where 

fi(G N t) - I exp {G N r)dT = [exp(G^t) - I]/G N (7.8o) 

Jo 

f 2 {Gf}t)— f exp(GNT)(t — r)dr = [exp(Gst) — G N t — I]/Glj. (7.86) 
Jo 

In order to implement our algorithm we have to get an approximation of the domain D 
which includes the eigenvalues of G^. One way to get it is by doing a Fourier analysis of 
the constant coefficients operator. 

~ (7 - 9o) 

where 

a = max |a(x)|: 6= max |6(x)|; 

0<x<2ir V n 0<x<2 *■' V 

Let Vk be an eigenvector of Gjv; then 

Mi = (7-10) 

and 

Ajt = -^r(cos(A;Ax) - 1) - *'-^-sin(A;Ax) — c (7-11) 

Ax^ Ax 


c = 


max 

0<x<2t r 


K*)l; 


or c = mm 

0<i<2»r 


|c(x)| 


(7.96) 
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is an eigenvalue of Gn • 

The tables in this subsection present numerical results for the fully descrete solution 
of (7.1) where 


k = 3 



a(x) = 1-/(2 + cosx) 

(7.12) 


b(x) = l./(2 + sin(x)) 

(7.13) 


c(x) = — 20./(2 + cosx). 

(7.14) 

From (7.11) we get 

A < Re Xk < B 

(7.15) 


\ImX k \ < C 

(7.16) 

while 



A= - 

~ Ax 2 20i B= 20/3 ; C= Ax' 

(7.17) 

Since 

\B-A\ » C 

(7.18) 

taking D as 

D = { x\A < x < B} 

(7.19) 

will be a relatively good approximation to the domain D. We have 



p(b) = (B — A)/ 4 . 

(7.20) 


Hence, in order to implement the new algorithm we have to interpolate the functions 

fi{pzt) = [exp (pzt) - 1 \/pz 
f 2 {pzt) = [exp (pzt) - pzt - 1 }/{pz) 2 

at points on ^D. 

In this case we can take the M interpolating points as the extreme of the scaled and 
translated Chebychev polynomial of degree M [Rivl74[. Thus 

% = \\{B - A)ij + B + A] 


j — 1 , . . • , M 
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where 


Xj = cos 


U - 1)tt 


M- 1 

In order to avoid roundoff errors one should arrange the interpolating points as described 
in Appendix A. Hence 


z i = j I K 5 ” A ) x > + B + A \ 


j = 1,... ,M 


where 


x 1 = l 


x 2 = -1 


(7.21) 

(7.22a) 


Xj = Re(u> 2 j- 3 ) j = 3, . . . ,M 


(7.226) 


and 3 are as in Figure (A.l) for the case M = 7. 

Since zy are real we can use algorithm 1 (3.2). The next two tables present the 
numerical results using algorithm 1 with zy defined by (7.21)-(7.22). 

Table 1. 

Mesh refinement chart - problem (7.1) 

Using algorithm 1. t = 1 


N 

M 

L 2 In 

32 

24 

.1108-01 

64 

52 

.2592-02 

128 

112 

.7033-03 


N - Number of grid points. 

M - Number of matrix-vector multiplications (Each evolution operator is approximated 
by a polynomial of degree M/2) 

L 2 I n - L 2 Error at t = 1. 

For sake of comparison we present in Table 2 a similar chart while using a standard 
second order in time scheme 

A 

U n+1 =u n + A tG N u n + -r-G%u n 

Jj 


(t = n • At). 


(7.23) 
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Table 2. 

Mesh refinement chart - problem (7.1) 
Using (7.23) t = 1 


N 

M 

L 2 T a 

32 

102 

.1108-01 

64 

408 

.2729-02 

128 

1632 

.6828-03 


(In this case M/2 is the number of time steps.) 

Comparing Tables 1 and 2 we observe that while in the standard second order scheme 
(7.23) M is proportional to N 2 (which results from the fact that At = 0(Ax 2 )) for the new 
algorithm, M is “almost” proportional to N. This phenomenon is explained and proved 
in [Tale 85]. For t — 20 we have the following results; 

Table 3. 

Mesh refinement chart - problem (7.1) 

Using algorithm 1 t — 20 


N 

M 

L 2 I n 

32 

24 

.1322-01 

64 

52 

.3407-02 

128 

112 

.1180-02 


Observing Tables 1 and 3 we notice that M does not depend on t. It can be explained 
as follows: For large t, || exp(Gjvi)|| is much smaller than the error which results from the 
space descretization. Thus, since (7.8) we have 

= —1/Gn 

f 2 (G N t) = -(G N t + I)/G 2 N 

which means that for large t, approximating /i(Gjvt) is equivalent to inverting G^,and 
approximating f(G N t) is equivalent to inverting G 2 N . 












27 


We did not implement the second order algorithm (7.27) for t = 20, but it is straightforward 
(by stability considerations) that for N = 32,64,128 we have M = 2040,8160,32640 
respectively. Hence for t — 20 and N = 128 the new algorithm is more efficient than the 
second order scheme by the impressive factor of 32640/112 = 290. 


II. Viscoelastic Wave Propogation. 

We have worked on this model with the Geophysical group in Tel- Aviv University 
headed by Dr. Dan Kosloff. A detailed report on the results for a general 2-D problems 
will be published elsewhere. In this subsection we describe the implementation of the 
algorithm for the simple 1-D model 


where 


U t = GU 0 < X < L 



° 

1 

0 \ 


/« 1 \ 

G = 

m u V 

0 

V 

; u = 

u 2 


a 

0 

— 1/t) 


\U 3 J 



m u 


= c 


1 - (1 ~ 



c — 2000 (the speed of sound) 


a = c 2 {l-r/rj)/rj 


(7.24) 


(7.25) 

(7.26 a) 

(7.26 b) 
(7.26c) 
(7.26d) 


u i is the pressure, U 2 is the pressure time derivative , uz is a memory variable and r, r\ are 
parameters of the problem. The initial data are 


ui(x, 0) = exp 




(7.27a) 


u 2 (a;,0) = « 3 (x, 0 ) = 0. 


(7.276) 


28 


The space descretization is done by Fourier method (GoOr81). Thus the semidescrete 
representation is 


(U N ) t = g n u n 

u n{x, 0) = Pjv[exp(— ^(x- ^L) 2 ] . 

Ulr{x,0) = Ufr{x,0)=0. 
Pn is the pseudospectral projection operator: 

N 

P Nf = a k exp[?(2 nk/L)x] 

k— — N 


(7.28 a) 
(7.28 b) 


(7.29) 


2N-1 


a k - 


2NC k 


^2 f{xj) exp[-i(2Trk/L)xj] 

3=0 


' G k = 1 for \k\ # N 

< 

^ Cn = C-n = 2 


(7.30) 


Xj = j Ax 

i = 0,...,2iV-l (7.31) 

Ax = L/2N 

and 

Gff = PfjGPff . (7.32) 

(Gjy- is a 6 N x 6 N matrix.) 

In our experiments we took 

N = 64 . (7.33) 

Ax = 20 . (7.34) 


Thus 

L = 2N x Ax = 2560 . (7.35) 

Since N = 64, Gpj is a 384 x 384 matrix. 

Approximating the domain D , which includes the eigenvalues of G^,i s done by making 
use of the following idea. Let us assume that the number of points K which are needed to 
resolve the coefficients of the operator G are relatively small. In this case, decreasing the 
space domain by a factor of K/N and using the same Ax gives us the matrix Gk- Since 
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K is small, one can use a library routine to compute the eigenvalues of Gk- The domain 
Dk , which includes this set of eigenvalues, is relatively a good approximation of D. In the 
present case is a constant coefficient matrix, and the domain D 4 achieved by computing 
the eigenvalues of G 4 (a 12 x 12 matrix) is exactly the same as D i28> The domain D 4 is 

D 4 = { z\z € [—a,0] or 2 G [— i&, if>]} (7.36) 

where a and b depend on r, 77. For this domain we have an analytic expression of the 
conformal mapping which maps the complement of the unit disc on the complement of D 4 
(5.12). The logarithmic capacity is given by (5.14). 

We ran two sets of numerical experiments. In the first one we took 

r = 0.1600890 x 10" 3 


and we have gotten 


r} = 0.1582262 x 10 -3 


a = 6320 ; b = 317 . 


Using algorithm 2 (3.26) for computing the solution at time level t = 0.1 results in 


N 

M 

L 2 I n 

128 

130 

.1588-02 

128 

150 

.2904-05 


(Since we do not have an analytic expression for the exact solution, we computed 
L 2 I n by comparing the numerical solution to another numerical solution achieved by using 
algorithm 2 with M = 300.) Similarly, using the second order in time algorithm (7.23) 
results in 


N 

— 

M 

L 2 T a [ 

128 

/?09 

.1902-1 


(For M < 632 the scheme is unstable.) 
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In the next experiment we took 

t = 0.1600890 x 10" 4 


rj = 0.1582262 X 10~ 4 


and we have gotten 


a = 63201 

The results at t = .1 were 


b = 317 . 


N 

M 

L 2 J n 

* 128 

400 

.2200-02 

128 

450 

.1096-04 


While using (7.23) results in 


N 

M 

L 2 T a 

128 

6320 

.1892-3 


(For M < 6320 the scheme is unstable.) 


III. Linear Systems. 

For the numerical experiments reported in this subsection, we have used a test matrix 
A NxN which is block diagonal. Each block is of the following shape 

U 2) • (7 - 37) 

Hence, the eigenvalues are 


A j = a,j it i\bjCj 


i=l,...,W/2. 


(7.38) 
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Two kinds of experiments have been carried out: 
bj - Cj => the matrix is normal. 
bj ^ Cj => the matrix is not normal. 

We solve 

Av = w (A is a 1000 x 1000 matrix) 

where 

v = [ui,...,t»iooo] T 

and 

v } r = (-1)*’ j= 1,...,1000 

and the initial guess v° is such that 


v°j = 0 j = 1, . . . , 1000 . 


Three methods have been tested: 


(1) I n - the method described in this paper 

(2) Cb - Chebyshev approach (Mant78) 

(3) M r - Minimum residual. 


The minimum residual algorithm is defined as follows: 
Given initial guess x° for A: = 0, 1, ... , until satisfied do 

r k = Ax k — b 

ajk = (r k , Ar k ) / (Ar k , Ar k ) 


(7.39) 


(7.40) 


(7.41) 
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1. First case - Cross shaped domain. 

In [SmSa85], the authors treat an example given by Hageman [HaYo81] which orig- 
inates from the solution of the neutron diffusion equation, where the eigenvalues of the 
Jacobi iteration matrix may be shown to lie on a the following Cross-shaped domain D: 



Re 2 


The conformal mapping from the exterior of the unit disc to the exterior of D is 

2 = 4>(w) = b + ^=^tx;2 + -L . (7.42) 

Since (7.42) we get that 


P = l/v/2 . (7.43) 

Thus the mapping function from the exterior of disc of radius p to the exterior of D is 

2 = xl>{w) =b+ + . (7.44) 

Since 

R = |0 _1 (O)| (7.45) 

we have 


R = 


b 2 + (b 4 — l) 1 / 2 


1/2 


Using (7.43) and (7.46) we get that the asymptotic rate of convergence is 

1 - 1/2 


p/R = [6 2 + (b 4 - 1 y/ 2 


(7.46) 


(7.47) 


According to the theory (4.1) 


I /(*) - Pm(z ) I < C(p/R) M ■ 


(7.48) 
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thus we can predict that after M iterations, the accuracy (in the normal case) will be 

L 2 Error(/ n ) £ ( p/R) M . (7.49) 

An improved prediction, which includes the constant C is: 

M 

L 2 Error(/„) = C{p/R)» Si [](1 - ±) (7.50) 

» = 1 *»' 

where z is any point at the domain D and z t - are the interpolating points. 

(In the next two tables we have chosen z = (6, 0) .) 

The numerical results are 


Table 7. 

Solution of ( 7.89 ) with b = 1.004 


M 

w«) m 

C{p/R)% 

Z/2Error ( I n ) 

L 2 Error (Cj,) 

L 2 Error (M r ) 




Normal 

N. Normal 

Normal 

N. Normal 

Normal 

N. Normal 

m 

1.575-3 

3.16-3 

3.87-3 

1.02+0 

1.85-1 

9.32+0 

9.67-2 

6.55+0 


2.8-6 

5.63-6 

6.9-6 

4.41-3 

1.04-1 

10.31+0 

5.37-2 

7.34+0 

H 

9.0-12 

1.8-11 

2.2-11 

2.29-8 

3.9-2 

7.65+0 

1.9-2 

7.13+0 
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Table 8. 

Solution of (7.S9) with b = 1.1 


M 

(p/R) m 

C(p/R)% 

L 2 Error (7 n ) 

L 2 Error ( Cb ) 

Z^Error (Af r ) 




Normal 

N. Normal 

Normal 

N. Normal 

Normal 

N. Normal 

10 

4.13-2 

8.99-2 

1.10-1 

9.69-1 

1.9-1 

9.9-1 

1.57-1 

8.0-1 

22 

9.03-4 

1.8-3 

2.22-3 

6.37-2 

5.01-2 

5.56-1 

4.07-2 

4.29-1 

42 

1.5-6 

3.1-6 

3.8-6 | 

1 

2.82-4 

6.34-3 

1.33-1 

5.1-3 

9.98-2 

82 

4.5-12 

9.0-12 

1.1-11 

2.15-9 

1.2-4 

4.81-3 

9.4-5 

3.27-3 


2. Second Case - Boomerang shape domain. 

Another example reported in [SmSa85] is Van der Vorst’s example. Let M be the 
matrix arising from the descretization of the P.D.E. operator u xx + u yy + /?ju x + 02 u y and 
let K be the incomplete CholesJaey factorization’, then the eigenvalues of the preconditioned 
matrix K~ l M, are sometimes observed to form the following boomerang-shaped profile 



Since, in this case we do not have an explicit expression for the conformal mapping, 
we have used numerical algorithm [TrefSO] for computing the interpolating points Z{. The 

predicted accuracy is 

M 

C(p/R)% = Ilf 1 - f ) 

,=i Zi 


(7.51) 
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while 

2 = (6,0). (7.52) 

Table 9 summarizes the results of this case. 


Table 9. 


M 

<+>/*)" 

L2ETTOT (I n ) 

i^Error ( Cb ) 

Z^Error ( M r ) 



Normal 

N. Normal 

Normal 

N. Normal 

Normal 

N. Normal 

20 

1.23-4 

5.4-4 

8.8-4 

1.18-1 

1.33-1 

2.42-1 

2.35+0 

40 

7.58-8 

2.36-7 

5.14-7 

4.3-2 

5.8-2 

1.13-1 

2.35+0 

60 

3.78-11 

1.03-10 

3.05-10 

1.71-2 

4.24-2 

5.59-2 

2.35+0 


In the next experiment we have shifted the domain to the left. 



According to the theory, the two methods: C<, and M r would not converge in this case. 
Table 10 verified this fact. 
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Table 10. 


M 


L2 Error (/ n ) 

I, 2 Error (Cj,) 

L 2 Error (M r ) 



Normal 

N. Normal 

Normal 

N. Normal 

Normal 

N. Normal 

20 

1.56-3 

6.65-3 

1.08-2 

1.15+0 

1.31+0 

7.09-1 

3.62+0 

40 

1.22-5 

3.8-5 

8.25-5 

4.93+0 

6.67+0 

6.76-1 

3.62+0 

60 

7.91-8 

2.16-7 

6.38-7 

23.0+0 

56.9+0 

6.58-1 

3.62+0 


3. Third Case - disjoint intervals. 

As mentioned previously (Section 2), since the complement of D is not simply con- 
nected anymore one should use a slightly different theory. In a future paper we will carry 
out a detailed analysis. 

In this subsection we report on a few experiments where D is 

D = hUl 2 (7.53) 

h = [ai,o 2 ] ; h = [o£ 3 , 0 : 4 ] (7.54) 

oci < a 2 < oca < 0 C 4 . (7.55) 


! 
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The error polynomial P N ( z ) is 

Pn (. 2 ) = Pn x ( 2 ) • Piv 2 ( 2 ) (N = N\ + N 2 ) (7.59) 

while 

(^) = n ^ 1 ~ z i z »■) ( 7 * 6 °) 

t=l 

N2 

p N 2 {z) = n^ 1 _ 2 / 2 ») • ( 7 - ei ) 

*=1 

We have 

l-Pv, (2)| - (pi/Bi) Wl 2 € lx . (7.62) 

Since the conformal mapping from the exterior of a unit disc to the exterior of an interval 


[a, 6] is 

x (u;) = 1 

^(5 - a)(u> + — ) + 6 + a 

2 to 

, (7.63) 

we get 

- 

I 

-O 

II 

a 

(7.64) 


R = 

= (a + 6+ \/a6)/4 . 

(7.65) 

Thus 


2 G 7i . (7.66) 

Pn 2 {z) satisfies 





|JV,(2)|<1 2 G /1 . 

(7.67) 

Therefore if 

/'v/cq-v/STN'' 1 _ c 

Vv/^+v«r/ 

(7.68) 

then 





|Pjv( 2)|<£- 2 G 7l . 

(7.69) 
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Now, for z £ I 2 we have 



max | P Nl ( 2 :) | ~ ka ” 1 1 = 1/ z \ j 

(7.70) 


\PsM\~W**) N ' 

(7.71) 

where 

. /7? - 

p 2 k 2 — . 

V®4 - v«3 

(7.72) 

Therefore, if 

kct4 l { P2 /R 2 ) N 2 <e 

(7.73) 

then 

|Pjv(*)| < e ze I 2 . 

(7.74) 


Using (7.66), (7.68), and (7. 73), we get that N 2 has to satisfy 

Iog(pi/i?ia 4 ) N log(l/fc) 

2 \og(p 2 /R 2 ) 1 \og(p 2 /R 2 ) 


(7.75) 


In Table ( 11 ) we report on experiments where cti, a 2 are fixed, and we increase 0 : 3 , <24 
such that (X 4 - <23 is constant. According to (7.75), it is easily verified that in this case 
N 2 — *• Ni . The predicted error is 


E r = max |Ptf (*0| = Pn{oc 4 ) = JJU - a *l z i) EE 1 “ 0:4 / 2 » ? ) • 

t=i »=i 


(7.76) 
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Table 11. 

ai = 1 <*2=3 tt 4 = a 3 + 2 



In the next set of experiments, we also increase the distance between 0:3 and 0:4. The 
results are reported in Table 12. 


Table 12. 

ai = 1 ai=3 0:4 = 3a 3 


<*3 


n 2 

E r 

L 2 Error (/„) 

L 2 Error (C&) 

L 2 Error (M r ) 





Normal 

N. Normal 

Normal 

N. Normal 

Normal 

N. Normal 

10 

10 

32 

1.6-7 


1.3-5 

7.7-7 

3.6-5 

3.8-4 

2.5-4 

20 

10 

38 

8.4-7 


9.7-6 

2.0-5 

1.1-3 

4.2-4 

1.5-3 

40 

10 

44 

2.2-6 

6.0-7 

9.8-6 

3.6-4 

2.2-2 

4.9-3 

1.3-2 

80 

10 

50 

3.9-6 

9.6-7 

1.5-5 

4.2-3 

2.8-1 

1.9-2 

2.5-2 


In the last set of experiments, we have negatives eigenvalues as well. In this case, one 
cannot use the two methods: Cb, M r . 











































40 


Table 13. 


ot.\ = 2 a<i 0:3 = — Q!2 0:4 = —ai 


<*1 

Ni 

n 2 

E r 

L 2 Error (J n ) 





Normal 

N. Normal 

2 

16 

16 

9.7-7 

3.5-7 

2.5-5 

4 

16 

16 

9.7-7 

3.5-7 

2.5-5 

8 

16 

16 

9.7-7 

3.5-7 

2.5-5 


8. Conclusions 

It has been shown that having an a priori knowledge on the distribution of the eigen- 
values of a matrix A, it is possible to construct an efficient algorithm for approximating 
f(A). The more accurate we locate the domain of the e.v., the more efficient is the al- 
gorithm. Two methods addressing this problem were described in Section 7. It seems 
to us that once this problem of finding the domain of the e.v. is solved satisfactorily, 
the algorithms described in the paper can be used as a numerical tool for many practical 
problems. 
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Appendix A. 

The set of interpolating points satisfy 

Zj = 4>{uj) j = (A.l) 

while w j are roots of the equation 

w N = 1 . (A. 2) 

The problem is: how to arrange t Oj such that the roundoff errors will be minimum. 

In a future paper we will carry out a detailed analysis addressing this problem. Ac- 
cording to the solution described there^we have 

w 1 ,w 2 ,...,wn (A.3) 


equally distributed on the unit circle, where 


wi = 1 ; u>2 = — 1 


and each “new” point twy is chosen such that the partial set 

wi,w 2 ,...,Wj (A. 4) 

will be as equally distributed as possible. For example, when N = 12 we have 



In the case of algorithm (3.27), (A.3) has to be modified as follows 


(A.5) 


xvi,w 2 ,xv 3 ,w 4 ,. ..,xv 2 j-i,xv 2 j, ...,rv N 
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while 



u>i = i ; 

w 2 = -1 

(A.6) 

W 2 j 

= w 2 j-\ 

2 < / < — 
~ 2 

(A.7) 

and each “new” pair of points is 

chosen such that the partial set 




.,W 2 j-l,W 2 j 

(A.8) 


will be as equally distributed as possible. Thus, for the case N = 12 we have 
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